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ABSTRACT 

Collisions resulting in fragmentation are important in shaping the mass spectrum 
of minor bodies in the asteroid belt, the Kuiper belt, and debris disks. Models of 
fragmentation cascades typically find that in steady-state, the solution for the particle 
mass distribution is a power law in the mass. However, previous studies have typically 
assumed that the mass of the largest fragment produced in a collision with just enough 
energy to shatter the target and disperse half its mass to infinity is directly proportional 
to the target mass. We show that if this assumption is not satisfied, then the power law 
solution for the steady-state particle mass distribution is modified by a multiplicative 
factor, which is a slowly varying function of the mass. We derive analytic solutions 
for this correction factor and confirm our results numerically. We find that this cor- 
rection factor proves important when extrapolating over many orders of magnitude in 
mass, such as when inferring the number of large objects in a system based on infrared 
observations. In the course of our work, we have also discovered an unrelated type of 
non-power law behavior: waves can persist in the mass distribution of objects even in 
the absence of upper or lower cutoffs to the mass distribution or breaks in the strength 
law. 

Subject headings: Asteroids - Collisional physics - Debris disks 



1. Introduction. 



Collisional evolution of many-body astrophysical systems in which the relative velocity between 
colliding objects is large compared to the escape speed and coagulation is unimportant is doni i nated 



by fragmentation. Examples of such systems include the asteroid belt (IDurda fc Dermott 



around young stars (jKennedy &: Wyattl . 



2011 



Bottke et aUboOfj ^. the Kuiper Belt dPavis fc Farinellal. Il997l: IPan fc Sari [20(151^. and debris disks 



Kenyon fc Bromleyl . l201ol ). As the large bodies are 



1997 



slowly ground down, a collisional cascade is launched, in which mass flows unidirectionally to 
smaller objects until at some scale it is flushed out of the system by some removal process (e.g. by 
Poynting- Robertson drag, radiation pressure, or gas drag). 
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In real astrophysical systems there is normally a large dynamic range between the scale TfiiYij 
at which mass is injected into the cascade and the scale rrirm at which mass is removed from the 
system. The mass rriinj can be defined as the mass for which the collisional destruction timescale is 
comparable to the age of the system. In highly evolved syste ms, such as the Kuiper Be l t, thi s scale 
corresponds to the characteristic mass of the largest bodies (jPan fc Saril . l2005l : iFraserl . l2009l ) . The 
collision time for the smallest bodies is usually much shorter than for the largest ones, so a steady- 
state can be set up for m^m *C m, ^ Jrimj. In such a steady-state, the mass distribution of particles 
evolves on the collision timescale of bodies at the injection mass scale and can be considered static 
on shorter timescales. 



Dohnanvil (|l969l ) was the first to construct a model of a fragmentation cascade that aimed 
at explaining the mass distribution of objects in the asteroid belt. He assumed that the internal 
strength of colliding objects is independent of mass with the implication that mBimt), which we 
define as the mass of the smallest projectile capable of dispersing one half the mass of a target 
of mass nit to infinity, has its mass linearly proportional to mt- Dohnanyi also assumed that 
mQ{mt,mp), the mass of the largest fragment produced in a collision between a target of mass 
mt and a projectile of mass m^, scales linearly with nit and is independent of rup. Under these 
assumptions, he showed that a fragmentation cascade allows a steady-state power law solution 

n{ni) (X ni~°^ , a = —11/6 (1) 

where n{m)dni is the number of objects in the size distribution with mass between ni and m + dm. 



Later, iTanaka et al.l (jl996l ) generalized Dohnanyi's result by assuming a self-similar model of 
fragmentation, which again entails msimt) oc mt, but now mQ(mt,mp) oc mtq{mp/mt), where q is 
an arbitrary function. They confirmed the power law form of the mass spectrum and showed that 
the value of a is determined by the mass dependence of the collision rate and that a reduces to 
11/6 if the collision rate is proportional to m/ (geometrical cross-section with mass- independent 
relative velocities). 

In reality, however, fragmentation does not have to be self-similar because the minimum energy 
necessary for disruption is not always linearly proportional to the target mass. For instance, if an 
object's internal strength is dominated by gravity, then Q^, the energy per unit mass required to 



5/3 



which is the 



shatter an object and disperse half o f its mass to infinity, scales as mtQ*p oc m^ 
case for objects larger than ~ 1 km (IBenz Asphaugl . Il999l : iHolsappld . Il993l : iBenz et ai 1 1994 ) 



O'Brien fc Greenberei (j2003l ) considered a model in which scales as a power law in m t , but took 
mo linearly proportional to mt- Under these assumptions, O'Brien Greenberg ( 20031 ) found that 
a steady-state power law solution for n{m) still exists, but that a differs from 11/6 unless is 
constant, even when the collision rate scales as mf^^. 



The steady-state power law solutions of lDohnanvil (j 19691 ) , iTanaka et al.l (| 19961 ) , and lO'Brien &: Greenberg 



(j2003l ) have since been confirmed numerous times by simulations, which have also shed light on 
non-power law effects present in astrophysical fragmentation cascades. These effects typically 
manifest themselves as waves superimposed on top of the steady-state power law solution and 
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are caused either by non-collisional ma ss sinks, e.g. due to the removal of small particles (S, 1 



p,m for L = L(^) by radi ation pressure (jXhebault &: Augereaul . 120071 : ICampo Bagatin et al.1 . Il994 



(O'Brien &: Greenbers. 


2003, 


2005): or a 


3 distribution ( 


Eraser . 


2009; 


Pan &; Sari, 



2005 



Kenyon fc Bromleyl . l2004l ). 



In this work, we describe a new source of non-po wer law behavior in fragmen tation cas- 
cades. We consider a model similar to the one used by lO'Brien fc Greenberei (120031). bu t with 
a more general form for mp. More specifically , previous researchers (IPetit &: Farinella 



O'Brien k Greenberd . l2005l : IWiUiams k Wetherilill994l : lKobavashi k Tanakal . l201C : 



1993 



de EKa k Brunini 



20071 ) have typically assumed that niQ = mtp{Ecoii/mtQ*jj), where EcoU is the kinetic energy of the 



colliding particles in the center of mass frame, and p is a function that varies from author to author. 
We instead consider the more general dependence tuq = m^p{Ecoii/'mtQ}j), which is motivated in 
f|2j We find that unless /i = 1, there is no steady-state power law solution. Instead, n(m ) is de - 
scribed by a power law with the same power law index as found by lO'Brien k Greenbergj (j2003l ). 
but multiplied by a slowly varying function of mass (i.e. n(m) oc m~"(p{m), \dln(p/dlnm\ <C 1). 
The non-power law effects caused by the slowly varying function show up as a smooth deviation 
from power-law behavior, which is quite different from the wave-like features described earlier. This 
deviation is significant when extrapolating over many orders of magnitude in mass, as demonstrated 
in D 



In the course of our investigation, we have also discovered that it is possible for waves to 
appear and persist in a collisional cascade, even if the particle size distribution does not contain 
an upper or lower mass cutoff, and the strength is described by a pure power law with no breaks. 
This is a completely independent type of non-power law behavior from the one caused by mo not 
proportional to m^. In astrophysic al systems, such waves ma,y be triggered by stochastic collisio n 



events between large planetesimals ()Kenvon k Bromlevl . l2005l : IWvatt k Dentl . l2002l : IWvattl . l2008l ) . 
However, the main focus of our paper is on the non-power law behavior that results when mo is 
not proportional to nit, and we defer a detailed discussion of these waves for the future. 

The paper is organized as follows. In ^ we introduce the general equations describing frag- 
mentation and discuss specific assumptions relevant to our model. In ?}3]we demonstrate that pure 
power law solu tions for a fragmentat i on ca scade are indeed possible if mo oc mtp{Ecoii/mtQ'^), 
consistent with iKobavashi k Tanakal (j2010l ) . We show in f^] that if the assumption mo oc rrit is 
not satisfied, then solutions are given by the product of a power law with the same index as in 
the mo oc nit case, and a slowly varying function of mass. We then find analytic solutions for this 
slowly varying function for monodisperse (all fragments having the same mass) and power law frag- 
ment mass distributions. We find that in the monodisperse case, the solutions for the steady-state 
distribution are not unique and can support waves. We confirm our analytical results numerically 
in ^ and discuss their validity and possible applications in ^ 
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2. Basic setup 

The number den sity distribution of particles in mass space, n{m), obeys the continuity equation 



(|Tanaka et alJ . Il996l ): 



dn{m)m ^ dF{m) 



dt 



dm 



0, 



(2) 



where the mass flux, F{m) is defined to be the amount of mass which flows past a point m in mass 
space per unit time. We wifl consider the evolution of the mass spectrum for m <C minj in which 
case we can adopt the steady-state assumption. Under this assumption, F{m) is constant, so there 
is no accumulation of particles at any scale. 

In order to write down the explicit form of F{m), we need to make some deflnitions. Dis- 
ruption of a target in a collision produces a spectrum of fragments characterized by the func- 
tion g{mf,mt,mp), where g{mf,mt,mp)dmf is the number of fragments in the mass interval 
{mf,mf + dnif) coming only from the target in a collision between bodies with mass rup and 
mt- Mass conservation then requires that 



g(mf,mt,mp)mfdmf = rrit. 



(3) 



Kobayashi Sz Tanakal pOlOl ) have shown that erosive collisions provide the dominant contri- 



bution to the mass flux, and in order to take them into account, it is useful to split g into two 
components 



g{mf,mt,mp) = grem{mf,mt,mp) + gej{mf,mt,mp). 



(4) 



Here, gej is the contribution to the mass flux from the continuous distribution of fragments ejected 
from the target and is normalized such that 



oo 

/ 



gej{mf,mt,mp)mfdmf = mej{mt,mp) 



(5) 



where rriej is the total amount of mass ejected from the target and dispersed to infinity. In 
catastrophic collisions with nip ^ ttt-r, rupj = mt. Howeve r, in collisions with mp < ms, the 



core of the target is left almost intact (iFujiwara et al.l . 119771 ). leaving behind a remnant of mass 
iT^rem = "nit — mej ■ This yields 



9rem 

{mf,mt,mp) = 5{mf - mremimt,mp)). 



(6) 



When mp = ms, then by definition of m^, farem = ^ej = 'rnt/2. A collision is commonly referred 
to as catastrophic when m^j > mrem and as erosive when m^j < mrem- We also clarify that ttiq is 
the largest fragment in the continuous distribution of ejecta, gej, to which mrem does not belong. 
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We next define f{m,mt,mp) as the mass fraction of debris with mj < m, which comes from 
the target only in a collision between a target of mass nit and a projectile of mass mp: 



f{m,mt,mp) = -^j g{mf,mt,mp)mfdmf. 




(7) 



This is consistent with lTanaka et alJ (|l996l ) and iKobayashi &: Tanakal (|201Cl ) and allows for projec- 
tiles that are larger than targets. We next split / into two parts, just as we did with g: 



f{m, nit, mp) = fej{ni, nit, nip) + frem{m, nit, nip). 
Using the definitions ([5|), ([6|), and ([7|), we have 



(8) 



fej{ni,mt,nip) = ^ j gej{mf,nit,nip)mfdnif, 





and 



frem{m,nit,nip 



{mt,mp)/mt , m > nirem{mt,nip) 
, ni<nirem{'mt,nip) 



(9) 



(10) 



We can now write down the mass flux in the form 



F{m) 



dmtnitn{mt) J dnipf{ni,nit,nip)n{nip)TZ{nit,nip), 




(11) 



where TZ{nit, nip) is the collision rate between bodies with mass mt and nip. The utility of splitting 
/ into two components (Eq. [8]) will become apparent shortly in §2.21 



2.1. Some simplifications 

When gravitational focusing is unimportant, the collision rate is given by 



TZ{nit,nip) = 




where (v) is an averaged collision velocity, which can be a function of mt and nip. To limit the 
number of parameters in our study, we will assume 

(v) = const, (13) 

so the collision rate becomes 

n{mt,mp) a (m^^ + m^^) ^ , (14) 
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but our results are easily extended to the forms of TZ considered bv lTanaka et al.l (|l996l ). who varied 
the power law dependence of TZ on mt and nip. 

It is natural to expect that disruption of targets with mass mt is dominated by collisions with 
projectiles having masses near or below the breaking threshold mB{mt). This is because the cross- 
section for catastrophic colli sions (defined by nip > msimt) (©) is dominated by the smallest 



particles as long as a > 5/3 (IDohnanvil. 



for m,p <C niBinit) (jKobavashi &: Tanakal . |2010| ). We now assume 



19691) ■ and the mass flux from erosive collisions drops off 



niB{mt) < mt, (15) 

which is typically valid for astrophysical fragmentation cascades. Then, TZ{mt) oc m/ for collisions 
that are responsible for the majority of the mass flux. This allows us to rewrite (jlip in the following 
form: 



F{m) oc J dmtmtn{mt)TZ{mt) j dmpf{m,mt,mp)n{mp). 



(16) 



The normalization of F{m) is unimportant since the only thing that matters for a steady-state 
solution is that F{m) is constant. 



2.2. Fragmentation Model 



Experimental data on collisional breakup (iGault &: Wedek indl ll 969l:lFuiiwara et al.l Il977|) and 
numerical simulations of high- velocity collisions (jBenz et all . 1 19 94; B enz fc Asphaug . 1999) suggest 
that the mass spectrum of fragments ejected from the target in a single collision can be reasonably 
well fit within a broad range of masses by a power-law with a cutoff at mo(mj, rripjl]: 



gej{mf,mt,mp) oc 



m 



, mf<mo[mt,mp) 
, mj > moimt,mp) 



(17) 



Equation ()17p is difficult to analyze for arbitrary dependencies of r]{mt,mp) and mQ{mt,m 

^ 1__ : I:j2__i:_._ x:_„x_J cir, oiil„ _x \ J^as the fomj^ 



Thus, we make the simplification, motivated in ^2.31 that 



gej[mf,mt,mp 



mej{mt,mp) f mf Ecoii{mt,mp)\ / ^^ /1c^ 

gej{mf,mt,mp) = — —ip -, — — i- , mo^B{mt) = mo{mt,mB{mt)). [18) 

m^j^{mt) \mo^B[mt) rntQci^t) J 



In the limit mp <^ mt, 



Ecoii{mt,mp) « mp{vf/2. 



(19) 



Fuiiwara et aP l| 19771 ). iTakaei et al.l (flQsi l. lDavis fc RvanI l|l990h find that a two or three slope power law better 



fits the data. 

^In some sense, this is more general than the form (|17|) . because g^j does not have to be a power law. 
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which together with the assumptions (jlSp and the definitiorjfl 

Ecoiiimt,mB{mt)) = Q)^{mt)mt (20) 

impHes 

E^oll{mum,)_ m, ^^(^^) « (21) 



mtQ*jjimt) niBimt)' 
It now helps to define the variables 



X = m/mo^Bimt) (22) 
y = mp/niBimt) (23) 
z = m/rrit (24) 

Then, with the form of gej given by Eq. (jlSp . fej becomes 

fej{m,mt,mp) = — ^^{x,y). (25) 

mt 

The normahzation of ^ is such that (,{oo, y) = 1, for y > 0, which fohows from Eq. ([5]) and Eq. ([9]). 
The prefactor in Eq. (|25p can be written as 

(26) 

mt mt 



mej{mt,mp) _ ^ _ mrem{mt,mp 
m- 

If we now make the assumption that 

mrem{mt,mp) 

-x{y), (27) 



mt 

we can absorb the prefactor in Eq. (I25p and write 

fej{m,mt,mp) = fej{x,y). (28) 

One can consider more general prescriptions for mrem/^t than Eq. (I27p . but the latter suffices to 
illustrate the non-power law behavior. It also follows from Eq. (f27l) that /rem has the form 



/„™(..rt= f • (29) 

[0 , z< x[y) 

so that finally we arrive at 

f{m,mt,mp) = frem{z,y) + fej{x,y). (30) 



Some authors (e.g. O'Brien fc Greenberj (2005)) assume that the projectile absorbs half of the energy, so 



Ecoii{''nt,mB{mt)) — 2QJjmt, but it does not matter which definition is adopted for our purposes. 



- 8 - 



2.3. Nonlinear scaling of rnQ si'^^t) 



We now address the natural q uestion of whether on e should expec t a nonlinear sca ling of 
rriQ^Bi'mt) in practice? Experimental ( Fuiiwara et al. . 1977 ) and numerical ( Benz &: Asphau ^ .^999) 
studies of collisional fragmentation suggest that the mass of the largest fragment formed in a high- 
speed catastrophic collision decreases wi th increasing collision energy. In particular, we focus on 
the experiments of iFujiwara et alJ (119771 ) . who fired polycarbonate projectiles of a constant mass 
and kinetic energy {nip = .37 g Vp = 2.6 km s~^) into basalt targets with masses in the range 
22 g < mt < 2900 g. They fit their data in the catastrophic regime with the following relation: 



max[mo{mt\mp,Vp),mrem{'mt\mp,Vp)] oc mt 



Ecoii{mt\mp,Vp 
mt 



7 ?s 1.24. 



(31) 



The vertical bar denotes the fact that the experiments of iFuiiwara et al.l (|l977l ) were performed 
at constant pr o jectil e mass and velocity. Since Ecou{mt\mp,Vp) const for mp <C mt, what 



Fujiwara et al.l (jl977l ) have shown is that mo{mt\mp,Vp) oc ml~^^ . However, we now make the 



following extension to their results. We assume Eq. (I3ip to be valid as a function of mp as well, so 
we write 



max[mo{mt,mp\vp),mremiint,mp\vp)] oc mt 



Ecoiiimtjmplvpj 



mt 



(32) 



In the highly catastrophic fragmentation regime {Ecoii ^ ^tQr))^ expect all of the fragments 
to be a part of the continuous fragment distribution, so that there is no remnant mass remaining 
("T-rem = 0). This means we can simplify Eq. (|32|) to the form 



mo{mt,mp\vp) oc mt 



E, 



coii[ mt,mp\Vp) 
mt 



(33) 



Using mp = msimt) in Eq. (|19p and making the usual assumption (|13|) then yields 

mo^B{mt) _ f mBimt) \~^ 
mt \ mt J ' 

Thus, unless oc mt, m^^s is not proportional to mt. Moreover, we see that if ms varies as a 
power law, then mo,_B varies as a power law as well. One assum ption that we h a ve ma de in deriving 
Eq. is that Eq. ([33]) is valid for mo < mt, even though iFuiiwara et al.l (jl977l ) obtained the 
power law relationship (Eq. (13ip ) by fitting primarily to data in the regime mo <C mt- Nevertheless, 
our analysis has shown that a power law dependence for mo,_B is plausible, and we discuss the matter 
further in ^6.21 



We now deduce what we would expect for th e power law exponent of mn.Bf m/:) in the range of 
target masses cq r isider e d bvlFuiiwara et al.l (119771). Experime nts and simulations (jBenz k. Asphaug . 



199S 



Holsapple 



1993 



Housen et al 



1991 



Benz et all . 1 1994 ) show that Q*j-) is well-described over 



a large range of masses by the expression 



Qorrit 



s/3 



(35) 
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where the exponent s is different for the strength-dominated and gravity-dominated regimes. If (v) 
is constant as we have been assuming, then it foUows from Eq. ()2ip that 

mBimt) = Bm^, /3 = 1 + s/3, (36) 

and consequently 

"T.o,B("T't) = Cm^, /i = 1 - 7s/3. (37) 



From simulations of impacts into basalt with Vp = 3 km s ^, iBenz Sz Asphauei ( 1999) found that 



s = — .38 in the strength-dominated regime. Using 7 = 1.24, the value measured bv lFuiiwara et al 



(j 19771 ) yields = 1.16. This is a small deviation from mo,B oc mt (i.e. ^ = 1), but enough to cause 



noticeable effects for real systems as we demonstrate in ^6.3i 

We next motivate our form for g^j in ^2.2l bv assuming the more general form 

9ej[mf,mt,mp) — — ^ rri [ 7 ^,y ] l38j 

and showing that it reduces to Eq. (jlSp if mo{mt,mp) is given by Eq. ()33p (with the assumption 
Vp = const). Using Eq. ([19]) and Eq. (f3^ in Eq. (j33|) . we can write 

(TTl \ 
7^ • (39) 
mBimt)/ 

Substituting this expression into Eq. (j38p and using the definition of y from Eq. (j23p . yields 

Making the definition il){x,y) = y'^'^il)i{xy'^ ,y), we arrive at Eq. (jlSp . 

3. Povi^er Lavi^ Solutions 

We now look for a steady-state power law solution for the mass distribution. Plugging Eq. ([1]) 
into Eq. (jl6p . using TZ{mt) oc rn^^^ , and using the form of / given in Eq. (j30p yields 

00 00 

F{m) (X / dmtm^^^~"' / dmpm~" {frem{z,y) + fej{x,y)) ■ (41) 



We now demonstrate how the power law solutions previously derived in the literature follow 
from this equation and elucidate under what conditions they fail. 
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3.1. Dohnanyi (1969) and Tanaka et al. (1996) case 



Dohnanvil (jl969l ) and iTanaka et alj (|l996l ) assumed a scale-free model of fragmentation with 
f{m,mt,mp) = f{m/mt,mp/mt), which in Dohnanyi's case was stated simply as tuq = Cnrit with 
C constant, and niB = Bnit with B constant. From Eq. (j30p . we see that this is equivalent to 
assuming mo,_B ^ "t^t and is constant. Changing the variables of integration lo x = z = ra/mt 
and y = mp/mt in Eq. (j4ip . and using Eq. ([8]) we have 



F(m) oc m"/3-2a / dxx^^-^^/s / dyy"" f {x , y) . 



(42) 



Taking a = 11/6 results in F{m) being constant in agreement with lDohnanvil (|l969l ) and I Tanaka et al 
(|l996l ). and we discuss the conditions under which the integrals in Eq. (j42p converge in Appendix 
El We will subsequently call a fragmentation model with niQ B = Cmt and thb = Bmt a "Dohnanyi 
model" . 



3.2. O'Brien & Greenberg (2003) and Kobayashi & Tanaka (2010) case 



O'Brien &: GreenbergI ([2003) and iKobayashi &: Tanakal (120101 ) went one step further and con- 
sidered a power law dependence of the strength as given by Eq. (f35|) . If s = in Eq. (ISS]) (i.e. (3 = 1 
i n Eg . ()36p), then this reduces to the Dohnanyi model. At the same time, lO'Brien &: Greenberg 
(j2003l ) and iKobavashi &: Tanakal (|2010l ) still took mo,_B = Cnit, so in their case x = z = m/mt and 



mp/mB{mt). Changing variables again to x and y in Eq. (j4ip and using Eq. ([8]), we find 



F(m) oc m5/3+(i+«(i-") / dxx-^/^-d+^d-") dyy-'' f {x , y) . 



(43) 



The mass flux is independent of mass if 



a 



P + 8/3 
/3 + 1 ' 



(44) 



which was derived by lO'Brien &: GreenbergI (120031 ) and IKobayashi fc Tanakal ()2010l ). and the reader 
is again referred to Appendix Rl for t he conditions under which the integrals in Eq. (1431) converge. 
The a rgum ents of Pan fc Saril (2005 ) are analogous to the calculations of O'Brien Sz Greenberg 



(j2003l ) and IKobayashi Sz Tanakal (120101 ) , but their qualitative nature has rid them of the need to 
worry about the scaling of mo with m^. We will subsequently call a fragmentation model with 
"^■0,5 = Cmt and /3 7^ 1 an "OBG model". 
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3.3. Failure of the power law solution. 

We now demonstrate that the power law solution ([T|) does not in general make the mass flux 
completely independent of m for any a, unless itiq^b oc mt as in ^3.113.21 To make the calculations 
tractable, we assume a power law dependence for in the form given by Eq. (j36p . and for mo,s 
in the form given by Eq. (I37p . 



Using Eq. (|4ip . we make the definitions 

oo oo 

) OC J dmtrri'J^ " j drripmp °'fremiz, y) (45) 

m 

oo oo 

Fej{m) oc / dmtmY^~°' / dmpUip" fej{x,y), (46) 



where F = Frem + Fej. Because x ^ z, in contrast to the Dohnanyi and OBG models, we change 
variables to z and y for the remnant fiux and to x and y for the ejecta fiux. This yields 

1 oo 

FrUm) oc m5/3+(i+«(i--)|dzz-8/3-a+/5)(i-)y"dyy-"/,e™(z,y) (47) 



m/mo,s(m) oo 

Fe,(m) oc rn'^-Hs/s+Ci+^a-a)) f dxx-i-^"'(5/3+(i+;3)(i-a)) /"^yy-a^^.^^^^^ (48) 



As before, a is given by Eq. ()44p in order for the m dependence outside both of the integrals to 
vanish. Now, however, m appears in the upper limit of integration in the integral over x in the 
expression for Fej, and unless mo,_B(?7i) oc m (i.e. ^ = 1), F{m) is not independent of m. This is one 
of the key conclusions of this work, and in the rest of the paper we will investigate the non-power 
law behavior of fragmentation cascades in detail. 

To keep things simple, we will assume that t he mass flux from reninants is negligible so that 
F{m) = F(,j{m). This is consistent with Fig. 3 of lKobayashi &: Tanakal ([20101), which shows that 
Fej /Frem ~ 10 when Q*j-) is constant. 



4. Non-Power Law Behavior 

From Eq. (j48p we see that the mass flux corresponding to the power law solution of the OBG 
model depends only weakly (logarithmically) on m. This motivates us to look for a solution of Eq. 
(jlOp in the form 

n{m) oc m~"ip{m), a = ^ (49) 
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where (p{m) is a slowly varying function of m: 

dip 



dm 



(50) 



This property can be verified a posteriori, after the explicit form of (p{m) is obtained. 
With n(m) given by Eq. (I49p . Eq. (j4ip becomes 

oo oo 

F{m) (X dmtm^/^~°'ip{mt) dmpm~"ip{mp)fej{x,y). 



(51) 



As discussed in Appendix [Aj the value of y~'^ fejix,y) typically drops off below y ^ k and above 
y ^ k, which means that m~'^ fej{x,y) is peaked at nip ~ kniBimt), where /c is a constant. In the 
case when erosion is neglected A; ~ 1, since collisions with projectile s of mass nip < uib co r itribu te 
no mass flux, but if erosion is included, then A; ^ 1 (Fig. 6 of iKobayashi &: Tanakal (120101 )). 
Together with the condition that (p(ni) is a slowly varying function of m (Eq. (I50p ). this allows us 
to expand ip{nip) in a Taylor series about nip = kmsinit): 



Lp{nip) K, ip{kniB{nit)) + 



dip 



din nir 



In 



With this approximation, the inner integral of Eq. (j5ip becomes 

oo 

d\n.ip 



ip{kniB{nit)) I dnipnip°'ip{nip)fejix,y) 



1 + 



m,p=kmB{mt) \kniB{mt) 



In 



(52) 



din nir, 



mp=kmB{mt) \kniB{nit) 



nir, 



(53) 



From Eq. ()50p . \d\nip/dlnm\ <ti 1, so as long as the peak in nip" fej{x,y) at nip = kniB{nit) 
is sharp enough that ln{nip/kmB{nit)) ~ 1 over the width of the peak, then the second term is 
negligible in comparison with the first. Thus, up to terms of order dlmp/dlnm, Eq. (j5ip becomes 



oo oo 
F{ni) <x J dnitni^^^~°'ip{nit)p>{kniB{nit)) j dnipni~'^fej{x,y). 



(54) 



Changing the inner variable of integration from nip to y, using the value of a from Eq. (|49p . and 
using the definition of niB{nit) from Eq. ()36p . we have 



F{ni) (X I dmtni^^ip{mt)ip{kmB{nit)) J dyy "fej{x,y). 





(55) 



Finally, defining 



fej{x) = / dyy '^fej{x,y), 
Jo 



(56) 
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and introducing the auxiliary function 

Q{mt) = ip{mt)(p{kmBimt)) . (57) 

we arrive at 

oo 

F(,„,«/.™„„r'e(™.)/.,(.) (58) 

m 

with X given by Eq. (j22p . Equation (j58p is the master equation for the two-step determination of 
ip{m): 

• First, given the exphcit form of fej{x) one needs to solve this integral equation under the 
assumption F =const to determine the behavior of the auxiliary function Q[m). 

• Second, having obtained Q{m) one must solve the functional equation, Eq. ([57|) . to determine 

We now perform this procedure explicitly for two specific forms of fej{x) — monodisperse ( ^4.ip 
and power law ( ^4.2p . 

4.1. Monodisperse Fragment Mass Distribution 

We first consider the special case of the monodisperse fragment mass distribution, which puts 
all fragments at a single mass scale mo{mt,mp) = mo^B{mt): 

gej{mf,mt) = — -6 {mf - mo,B(mt)) . (59) 

This singular fragmentation model can be thought of as a very crude qualitative approximation to 
any fragmentation law that has most of the debris mass concentrated at one scale. It allows us to 
obtain some interesting analytical results and serves as a simple stepping stone for the more general 
case considered in fHl 



The fragmentation law (|59p implies 



/e,(x) = r' ""^l, (60) 
[1, X> 1 

where x = m/mo,B(mj) ( ^2.2p . Plugging this into the master equation (Eq. ([58]) ) one obtains 

"io,s(»n) 

F(m)oc [ —e{mt), (61) 
J mt 

m 
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where rhQ^Bi'm) is a new function defined as a function inverse mo,_B(?7i) (i.e. ''7io._B(?7io,B(m.)) = m). 
Upon differentiation with respect to m, expression (j6ip results in 

e(m) e (mo Bim)) dmo Birn) 

= — ~ — r~\ A ' (62) 

m rriQ^Bvn) dm 

where we have used the fact that F{m) is constant. This functional equation is valid for arbitrary 
rriQ^Bi^t) provided that the fragmentation law is monodisperse. 

We now focus on mo,B(m) in the form 

mo,B{m) = ( — ) , (63) 



.CJ 

which is valid for uiq^b given by Eq. (j37p . Plugging this expression into (j62p we find 



e,™,4e((g)-). (.4, 



Introducing the new variable t = ln(m^ ^/C) = lii(m-/mo,B("i-)) and the new functioij^ 01 (t) 
G ({Ce^y^^^'^A, Eq. (IMl) becomes 



Qiit) = -91 (-) . (65) 

This has the form of a homogeneous functional equation 

q{au) = bq{u), (66) 
(o and b are constants) which has the solution 

q{u) = T {In u)u\ A = (67) 

ma 

Here T(s) = T(s + In a) is aii arbitrary periodic function with period In a, which can be constant 
(jPolvanin &: Manzhirovl . Il998l ) . This implies that the solution of Eq. (j65p is 

e,W = ^, (68) 

SO that finally 

T(ln(ln(m/mo,B("i))) , i ^ ran\ 

e(m) = — — , r(s) =r(s + ln^). (69) 

m(m/mo,B(mjj 

Clearly this solution is inapplicable to the OBG case of = 1, because the variable t then 
reduces to a constant. However, for = 1, Eq. (f64|l already has the form (|66|) . so that its solution 
is 

G(m) = r(lnm), r(s) = r(s + In C7). (70) 



^We are grateful to Jeremy Goodman for suggesting this transformation. 
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In particular, G(m)=const, and consequently 93(m)=const, is one of the possible solutions, so that 
n(m) simply proportional to m~" wit h a given by (I44D is a viabl e solution for a fragmentation 
cascade with fi = 1, in agreement with lO'Brien &: Greenbera (j2003l ). 



However, the existence of periodic solutions brings about the possibility of having waves in 
the mass distribution of objects while still having F{m) constant, even for /i = 1. The pres- 
ence of waves at masses m/rrirm '• ^ 1 in fragmentation cascade s having a lower ma s s cuto ff has 
been previously demonstrat e d bv ICampo Bagatin et alj (119941 ') . iDurda k. DermottI (jl997l ). and 
Thebault &: Augereaul (j2007l ): lO'Brien &: Greenbergj (j2003l ) found waves to appear whenever the 
scaling of specific energy necessary for disruption Q*j~) with object mass changed abr uptly (e . g. due 



to an object's self -gravi ty becoming more importa nt than its internal strength); and iFraserl (|2009l ). 
Pan &: Saril (j2005l ). and lKenyon &: Bromley! (|2004l ) have shown waves to be present at the transition 



from a collisionally evolved to a primordial size distribution (i.e. at m/m 



11 



The nature of the waves we have found is different from the ones discussed by previous authors, 
since they exist even when rrirm = 0, rriinj = oo, and is given by a pure power law without 



breaks (Eq. ([36]) 1. In astrophysic al systems, these kinds of waves could be trigger e d in stochasti c 
collisions of large planetesimals (jKenvon Bromlevl . l2005l : IWvatt Dentl . l2002l : IWvattl . l2008l ). 
Most of the mass in such a collision would be in particles of size mo(mt,mp), and if the density of 
particles with mass mo created in the collision is comparable to or exceeds the local disk density 
of such particles, a wave will be triggered. However, a proper treatment of these waves needs to 
account for a non-monodisperse fragment mass distribution, which could damp them, so we leave 
this subject for future work (Belyaev & Rafikov, in preparation). 



4.2. Povi^er JjSw Fragment Mass Distribution 

We now assume that the mass spectrum of fragments produced in a collision g(.j{mf,mt^mp) 
is a power law with an index —r] having a cutoff at a maximum fragment mass mQ{mt,mp) = 
mo,_B(m.(). This allows us to write 



fejix) 




(71) 



where as before x = m/mQ^B{fnt) ( ^2.21). This fragmentati o n model resembles reality, since ob- 



servational and experime ntal evidence (iGault k, Wedekindl. 
as numerical simulations ( Benz et al . 19941 : Benz &: Asphaug . 

Qejim, mt, mp) at small fragmen t masses (^2.2D. Various flavors of such a power l aw model have been 



196S 



Fujiwara et al.l . 119771 ) as well 



19991 ) suggest power law behavior of 



'il (Il9e 



adopt e d in theoretical studies bvlDohnanvil (119691). IWilliams fc Wetherilll (119941') .lO'Brien &: Greenbere 
\200i ). kenvon fc Bromlevl (|20ld ). lOavis k Farinellal (|l997l ). IXobavashi Tanakal (|20ld ). etc. 



- 16 - 



Plugging Eq. ([7T]) into the master equation (Eq. ([58]) ). we find 

rho,B{m) 



Q[mt). (72) 



F(m) oc [ '^Qimt) + [ 
J rnt J 



mo,B(m) 



Differentiating this expression with respect to m (keeping in mind that dF{m) / dm = 0), dividing 
the resultant expression by m^~'^, and differentiating again one finds 

dlnmo,B(m) f \\ r,t \ ^ dQ{m) 

— ^ e (mo,ij(m)) = e(m) - — . (73) 

ainm 2 — ij dlnm 



This equation reduces to Eq. ()62p if one takes the limit t] — > — c« which corresponds to all of the 
fragments' mass concentrated in objects with mass mo,s, and is thus equivalent to a monodisperse 
fragmentation law. 



Until now, our treatment was rather general and based solely on assumption (|7ip so that the 
functional equation (j73p is valid for any form of rriQ^B- We now take mo,B in the form of Eq. (j37p 
and find 

e(mo,B(m)) ^ ^^^^ _ 1 dQjm) ^^^^ 
^ 2 — r] dlnm 

In the limit rj — > — oo, corresponding to the monodisperse case, the second term on the right hand 
side of ([7i|) is zero, and the solution is easily verified to be 

®("^) = w / ^ / ^^ • (75) 
ln(?n,/mo,_B(?n.)) 



This solution could also have been obtained from Eq. (|69|) by setting T = 1. 

For the non-monodisperse case, we can solve Eq. (|74p by first introducing a new indepen- 
dent variable w = 1/ ln(m/mo,_B(m)), and the new function Qi{w) = Q{m) = Q({Ce^^'^)^^^^~'^')). 
This is a natural choice, since G(m) = w is the exact solution for the monodisperse case if 
T(ln(ln(m/mo,_B(?n.))) = 1. With these definitions we can rewrite Eq. (j74p as 

= id>i(w) w ; . (75) 

/i ^ ^ 2-r] dw ^ ^ 

We next look for the solution of this equation in the form of an infinite series 

oo 

Gi(«;) = J]W (77) 

k=l 

(note that = as long as fi ^ 1). By plugging this ansatz into Eq. (j76p and changing the index 
of summation in the last term we get 

oo oo ^ oo 

^ Akfs'-'w'' = Auw^ + ^ Y.^k - l)Ak.^w\ (78) 

k=l k=l ^ k=2 
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We can set Ai equal to any value and this corresponds to an overall normalization of 0i. We then 
obtain the other coefficients from the recursive relation 

Unfortunately, this series only converges for /i > 1. 

For fi < 1, we have found numerically in §5.1.11 that the formula 

ln(m/mo,B(m)) + 1/(2 - t]) 



gives good results up to < 1.7. Equation (|80|) is a version of the monodisperse solution ([75 
which has been shifted in Inm, and is accurate up to terms of O^w^) in Eq. ()78p . 



5. Numerical verification of non-power law behavior 

Having obtained solutions for G(m) for a couple of fragmentation models, we are now in a 
position to determine <p{m) from Eq. (j57p . The analytical calculations involved in this process 
are rather cumbersome and we refer the interested reader to Appendix [B] for the mathematical 
details. There, we describe the general method of solving for ip{m) which works for niB given 
by Eq. (I36p and for arbitrary 0(m). However, we provide explicit analytical results only for 
Q{m) corresponding to the monodisperse case in Appendix IB. 3 1 In this section, we compare these 
analytical results with numerical calculations of fragmentation cascades. The latter were carried 
out using a fragmentation code which is described in detail in Appendix O For simplicity, we 
ignore erosion in our calculations, which amounts to setting k = 1 m Eq. ()57|) . 



5.1. Results for /-i 7^ 1 and mB{m) = m. 

As discussed in ^2.1^ real astrophysical systems typically have rriBirnt) <^ nit. In order to 
qualitatively understand the non-power law behavior, however, it is instructive to consider a sim- 
plified model in which there is no erosion. Next, we also assume that a target can only be broken 
by projectiles that are of very nearly the same size as itself: 

mB{mt) = (1 - e)mt, e < 1. (81) 

Although such an assumption is unrealistic, it is useful for getting a qualitative picture of the 
non-power law behavior. 

In the model we have just introduced, which we will call the tub = ma model, the analog of 
Eq. dil]) is 

00 rhB{rat) 

F{m) (X J d'mtm^^^~'^ip{mt) j dmpmp'^ip{mp)fejix,y), (82) 

mB{irit) 
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where niBim) is defined to be the largest mass that can be broken by a projectile of mass m (i.e. 
mBiniBim)) = m). The power law slope of this model is a = 11/6 just like the Dohnanyi case, and 
it is straightforward to show that the master equation (Eq. ([58]) ) and Eq. ([57|) are both still valid 
for the niB = m model. Now, however, Eq. ()57p is trivial to solve, since we have mB{mt) ~ mt, 
which implies that Lp{m) ~ 0(m). Comparison between the numerical and analytical behavior 
of ip{m) then becomes a direct verification of our solutions for 0(m) obtained in Q and we avoid 
the intermediate steps involved in solving Eq. (j57p for arbitrary ms- 

In our numerical calculations we implement the case of mB(mt) = nit by taking mB{mt) = 
0.99mf. We specify the amount of time for which simulation was run in one of two ways, but both 
make use of the collision time defined generally as 

-1 

(83) 

Here, we have assumed a geometric cross-section for collisions and mass-independent relative ve- 
locities as in N2.1I First, we can specify how long a simulation was run in terms of the number of 
collision times of the smallest particles in the simulation. This is a natural unit of measure to use, 
since this is the timescale on which the low mass end of the particle mass distribution evolves, unless 
the initial distribution is already in a steady-state. Second, we can specify how long a simulation 
was run by the location of the collisional break, ^^)r'eafc5 

at the end of the simulation, t = tend, 
which is given implicitly by Tc{mbreak,tend) = tend- The quantity rrihreak is especially useful when 
interpreting simulation results for which the distribution was initialized to be in a steady-state to 
check an analytical steady-state solution ( ^5.1.2l and l5.2p . In this case, there should be no evolution 
of the distribution in time, and so no break should develop. However, rrihreak gives the mass below 
which we would have seen evolution, if the analytical steady-state solution were incorrect. 



Tc{mt,t) 



2/3 



fhB{mt) 



TTV 



dmpn{mp 



,t) (m] 



1/3 



+ m 



1/3 



5.1.1. Monodisperse Fragment Mass Distribution. 

We first consider the ansatz (18ip with p ^ 1 and the monodisperse fragment mass distribution. 
Using Eq. ([69]) we find 

(p{m) = ^ , (84) 

y^ln(m/mo,B(m)) 

where we have set the arbitrary periodic function T to a constant. 

To verify Eq. (j84p numerically, we initialize a pure power law mass distribution no{m) oc m~" 
with an index a = 11/6, which is the steady-state solution of the Dohnanyi model. We then expect 
that for fi ^ 1 the shape of n(m,t) will gradually evolve from ?io(m) towards the correct solution 
()49p with ip{m) given by Eq. (j84p . To illustrate this evolution in Fig. [T]and all subsequent figures, 
we plot the "effective population index" aeff{m,t) = —dlnn{m,t)/dlnm as a function of m at 
different moments of time. This way of representing the evolution of n(m, t) naturally highlights 
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Fig. 1. — Population index dlnn/dlnm = —Oeff of the particle mass distribution vs. mass for two 
different values of fi in the relation uiq^b = Cm'^. Both panels display the evolution of Og// for 
mB(m) = m starting from an initial distribution no{m) oc m~^^^^ (the Dohnanyi model). Panel 
(a) shows the smooth convergence of n{m,t) from nQ{m) to the solution n(m) = m~°'ip(m), with 
(p{m) given by Eq. ()84p and shown by the dashed line. The model parameters are mB{Tn) = m, 
fj. = 0.95, and mo,B(m- = 10^"*^^) = 0.3m (mass scale is arbitrary), which sets C in the expression 
for mo,B. The curves A-D show the solution after 25, 250, 2500, and 25000 collision times (defined 
in §5.ip of particles with mass m = 10~^*. Panel (b) shows the numerical solution starting again 
from a Dohnanyi distribution as in panel (a), but with = 0.7 and all other parameters the same. 
The simulation was run for 5000 collision times of particles with mass m = 10~^*. The numerical 
solution (solid line) no longer converges smoothly to the analytic solution (I84p (dashed line), but 
oscillates about it. Given the location of the leftmost peak (dashed arrow), Eq. (j69p correctly 
predicts the locations of the other peaks as indicated by the solid arrows. 



the non-power law behavior, since Oe// for n given by a power law as in Eq. ([T]) appears as the 
horizontal line Oe// = a in such plot. Thus, any deviation from a horizontal line is indicative of 
non-power law scaling of n{m). 

In Fig. [1^, we display the case of close to unity (^u = 0.95), which according to ([8^ 
corresponds to an almost constant since m/mo^Bi^) oc m*^"^; note the small range of variation of 
(Xeff in Fig. [T^. For this value of /i we indeed find that the initial power law distribution smoothly 
converges to the analytic solution (I84p over time. 

In Fig. [TJj, we show the evolution of Oe// for /U = 0.7, appreciably different from /x = 1. One 
can see that in this case the numerical solution does not converge towards the analytical solution 
in a smooth fashion. Instead the numerical solution oscillates in mass space about the analytical 
solution. The frequency of these oscillations is correctly predicted by Eq. ([69j) . and the positions of 
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the peaks computed according to this formula are shown by the arrows in Fig. [T]3. Their agreement 
with the locations of numerical peaks proves that waves can indeed be spontaneously generated 
and persist with no indication of damping in smooth fragmentation models (i.e. for functions niB 
and mo not having any breaks caused by the abrupt changes of the material properties of colliding 
objects), at least for a monodisperse fragment mass distribution. 

In this study, the appearance of waves is undesirable as it complicates the comparison between 
numerical and analytical solutions for n(m), and one would like to avoid their excitation. Comparing 
the two cases depicted in Fig. [1] suggests that waves get generated when the initial distribution 
no(m) is significantly different from the analytical, non-oscillatory steady-state solution. Motivated 
by this observation, we start from an initial distribution no(m) which is identical to the analytical 
steady-state solution in subsequent calculations. We then expect that the numerical solution, 
n{m,t), will not deviate from no(?n.) as time goes by if our steady-state solution is correct; if it is 
not, then n{m,t) will evolve away from nQ{m). 

Despite the complications related to the appearance of waves, it is clear that in the monodis- 
perse case, n(m) exhibits non-power law behavior for fi ^ 1, and our analytical solutions (j69p and 
()84p accurately describe the deviation from the pure power law solution. 



To better understand the qualitative behavior of the steady-state ae//(m) in Fig. [H we use 
Eq. ([37]), Eq. ([M]), and Eq. ([HI]) to write 



equal to the mass of the target (i.e. mo,_B(mo) = nT-o)- For < 1 a fragmentation cascade can only 
exist for m > tjiq, in which case a^ff > a, as can be seen in Fig. [T] (a = 11/6, since (3 = 1 for 
niB = fn). Thus, in the monodisperse case with ruB = rn, the /x < 1 collisional mass spectrum 
is steeper than in the Dohnanyi model. The deviation of Qe// from the Dohnanyi slope increases 
as m ^ niQ, and decreases for m ^ mg. For example, in the case shown in Fig. one has 
uiq = 1.8 X 10~^^ and Oe// deviates from 11/6 by ~ 0.12 already at m = 10~^^, i.e. at m/rriQ ~ 55. 

On the contrary, in the monodisperse case with ^ > 1, a fragmentation cascade is possible only 
for m < rn-Q, and Eq. ([85]) then predicts that a^ff < a. Thus, the collisional mass spectrum for 
/i > 1 has a shallower slope than the Dohnanyi solution. This can be seen in Fig. [5^ (curve labeled 
"— oo"). Fig. [3)3 (dotted curve), and Fig. (dotted curve). In that case, the biggest deviations of 
aeff from a are observed at large masses, as m — >■ ttiq. 



1 



(85) 




where mg = mass scale at which the mass of the largest fragment becomes formally 



5.1.2. Power Law Fragment Mass Distribution 



We now consider ansatz (|8ip with the power law fragment mass distribution explored in ^4.2[ 
In this case, ip{m) is given exactly as the square root of Eq. ([78p for > 1, and approximately 
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Fig. 2. — Population index of the mass distribution vs. mass for a power law fragment mass 
distribution and two different values of ^u. nQ(rn) was initialized to be the analytical steady-state 
solution for each model. The dot-dashed line indicates the a = 11/6 "Dohnanyi" power law index 
which would be expected for the = 1 case. The labels in each panel indicate the value of power 
law index r\ of the fragment mass distribution for different models, and a monodisperse model 
(corresponding to 77 = —00) is plotted for reference (its behavior is given by Eq. (j84p ). Panel (a): 
The numerical calculations were compared with the analytic series solution (Eq. (I78p ) truncated 
after 20 terms. Curves for two models both having 771^(771) = m, jj. = 1.5, mQ^B{fn = 10^^) = 
0.1m, and either 77 = 1.8 or 1.9 are shown. The analytic solutions (dashed lines) agree well with 
numerical runs (solid lines), although there is some deviation for the 7/ = 1.9 case. The location 
of ruhreak ( §5. ID in the simulations is mi^reak ~ 10'^. Panel (b): Similar to panel (a), but now the 
approximate solution (Eq. (j80|) ) was computed for three models, each having mB{m) = m, fi = 0.5, 
'^o,b("t- = 10^^^) = 0.1?Ti, and either 77 = 1, 1.7, or 1.9. There is good agreement between the 
analytic solutions (dashed line) and numerical results (solid lines) for tj < 1.7. Here rni,reak — 

10^ 

although we have truncated the mass range at tti = 1 to highlight the differences between the 
solutions. 



as the square root of Eq. (j80|) for /i < 1 . We again test these solutions numerically, but this time 
starting with the analytical steady-state mass distribution (as discussed in the previous section), 
since we are not interested in waves. 

For the fi > 1 case, we find that the solution given by the the series in Eq. (I78p converges 
quickly, and truncating the series after the first twenty terms gives a result which shows little sign 
of evolution for 77 = 1.8, (Fig. [2ti). Since in the 77 = 1.8 case numerical n{m,t) does not evolve 
significantly from the initial distribution given by the steady-state solution (Eq. ()78p ). this implies 
that our analytic solution for /x > 1 is indeed correct, even for 77 very close to 2. In Fig. [2]3 
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we show the analogous calculation for the /i < 1 case. We initialize no to be the approximate 
solution given by Eq. ()80p and find that for ij < 1.7 this solution works well. However, for larger 
values of r] deviations between the steady-state numerical and analytical solutions become apparent. 
Nevertheless, the general qualitative behavior is still reproduced by Eq. (jSOp even for r/ close to 2. 

It is easy to see from Fig. [2] that as rj ^ 2 the behavior of a^ff flattens out and approaches 
a power law solution with slope given by Eq. (j44p . which is just a = 11/6 for vtib = rn. Such 
behavior can be understood by looking at Eq. ()74p . which reduces to dQ{m) / dm = 0, in the limit 
of r/ — > 2. This means that 0(m) does not vary with mass in this limit, which implies that (/?(m) is 
also constant if mBim) = m. Thus, in the limit r] ^ 2 when fragments produced in an individual 
collision are uniformly distributed in Inm, the steady-state solution is a power law. 

On the other hand, the smaller r] becomes, the more the fragment mass distribution g is dom- 
inated (by mass) by the largest fragments, which makes it approach a monodisperse distribution. 
One then expects that when r] is reduced, Oe// should tend towards the monodisperse form given in 
Eq. (j85p . and this trend is clearly seen in Fig. [2l It is also worth noting that for a given value of fi, 
the largest deviation of cxeff from pure power law behavior occurs for the monodisperse fragment 
mass distribution, which is obtained in the limit r] —oo. 

5.2. Results for /U 7^ 1 and general mB(m). 

We now describe our results for the more realistic case when 7^ m^. uiq^b is again given 
by Eq. ()37p with fJ. ^ 1, which is necessary for the non-power law behavior. The general approach 
developed in Appendix iBl allows us to compute ip{m) for arbitrary 0(m), including that obtained in 
^4.2l for the power law fragment mass distribution. To keep things simple, however, we only provide 
a comparison between numerical and analytical results in the case of a monodisperse fragment mass 
distribution, for which we derived closed form analytical expressions in Appendix [Bl 

We start by looking at the case of /3 = 1, so we have mB{fn) = Bm, B < 1. For a given set 
of parameters, we initialize nQ{m) to be the analytical steady-state solution given by Eq. (IB23P if 
^ > 1 or by Eq. (1B25P if < 1. We then check whether n{m,t) evolves away from no{m). If it 
does not, then no(m) is the steady-state solution. 

In Fig. [3^1. we compare the analytical formula ()B25P for fi = 0.5 with the numerical solution. 
Despite small deviations between the two solutions (likely due to some problems with boundary 
conditions, see Appendix [Cl and, possibly, weak wave excitation) the overall agreement between 
them is quite good. Figure [Sb provides a comparison between our formula ()B23p for /i = 1.5 and 
the numerical solution. In this case, the two solutions agree with each other so well that they are 
hard to distinguish. 

We next look at the case of 7^ 1 and /3 7^ 1, so that ms is given by Eq. (p6|) . Proceeding 
as before, we display in Fig. |4] the evolution of numerical curves for Oeff away from analytical 
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mass mass 

(a) (b) 

Fig. 3. — Population index of the mass distribution vs. mass for mB{m) = Bm, B = 10~^, and two 
different values of /i. no(m) was initialized to be the analytical steady-state solution for each model. 
The dot-dashed line indicates the a = 11/6 "Dohnanyi" power law index which would be expected 
for the = 1 case. The dotted line represents the monodisperse solution with mB{m) = m. Panel 
(a): Analytic solution ()B25p . (dashed line) is compared with the numerical result (solid line) for a 
monodisperse fragmentation model having fi = 0.5, and niQ^Bi'm' = 10~^^) = 0.1m. The location 
of rrihreak IS at mtreak ~ 0.01. The effect of the imposed boundary conditions (Appendix [C|) is 
evident as the straight line between m = 10~^^ — 10~^®. Panel (b): Similar to panel (a), but now 
the analytic solution ()B23P was computed for a model with fj, = 1.5, and mo,B(?Ti = 10^^) = 0.1m. 
The dashed and solid lines are hard to distinguish here because the numerical solution does not 
significantly evolve away from the initial analytical solution. The value of m^reafc is n-breafc = 

solutions computed in Appendix lB.3.2[ The panels in this figure show our results for /3 < 1, /i > 1 
and /3 > 1, n < I. We do not show the results for /3 > 1, /u > 1 and /3 < 1, ^ < 1 due to numerical 
difficulties with imposing boundary conditions in these two cases (Appendix [C]) . The agreement 
between the analytical formula ()B27p for /3 < 1 and the numerical results displayed in Fig. is 
quite good, and the same is true regarding the agreement between the formula ()B26P for /3 > 1 and 
the numerical results displayed in Fig. Uja. 

In both Fig. [3] and Fig. HI we display Oe// computed for m^ = m and a monodisperse 
fragment size distribution ( §5.1.ip by a dotted line. One can see that even though we are now using 
m^ <C m, the solutions are qualitatively similar to the ms = m case. Thus, one can use the fully 
analytic solution (|i9]) with ip{m) given by Eq. ([8^ to get a quantitative picture of the non-power 
law behavior regardless of the precise form of m^- We also note that the solutions for Oe// shown 
in Fig. [3] and Fig. U lie above the solution corresponding to the ms = m case. Thus, niB ^ m 
gives rise to a function ip{m) which is shallower than for the uib = m case. 
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mass mass 



(a) (b) 

Fig. 4. — Population index of the mass distribution vs. mass for the cases /3 < 1, fi > 1 and 
/3 > 1, n < 1. Calculations were initialized at the analytical steady-state solution for each model. 
The dot-dashed line indicates the a given by Eq. (I44p — the OBG power law index which would 
be expected for the ^ = 1, /3 ^ 1 case. The dotted line represents the monodisperse solution 
with mB{m) = m. Panel (a): The analytic solution (IB27P (dashed line) is compared with the 
numerical result (solid line) for a monodisperse fragmentation model with (3 = 0.8, ^ = 1.5, 
nT'0,B{'m = 10^^) = 0.1m, and mB{m = 10~^^) = 10~^m. The location of the break is rrihreak = 10^. 
Panel (b): Similar to panel (a), but now Eq. ()B26P was used to initialize the analytic solution for 
P = 1.1, = 0.9, mQ^B{n^ = 10~^^) = 0.03m, and mB{m = 10^^) = 10~'^m. The location of the 
break is at miyreak = 10~^. In both panels (a) and (b), the numerical and analytical results are in 
agreement. 



6. Discussion 



6.1. Validity of |d In v9/d In m[ < 1 

In deriving our analytical results, we have assumed that \dlnip/dlnm\ <C 1 (Eq. [50|) . Results 
from the previous section demonstrate that the qualitative behavior of ip{m) is insensitive to the 
specific form of m^. Thus we can get a sense of when this assumption is valid by using the form 
of dlmp/dlnm for the monodisperse case with m^ = m: 

dlnm 21n(m/mQ)' 

In order to have \dlinp/dlnm\ <C 1, we must have | ln(m/mg)| ^ 1, where itiq was defined in ^5.1.11 
In practice, we find that even for | ln(m/mQ)| ~ 4 our analytical solutions give an accurate descrip- 
tion of the non-power law behavior. For instance, in Fig. [2]3 it is clear that for the monodisperse 
fragmentation law, the exact solution for ip{m) (Eq. ()84p ) works very well, even though we have 
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ln(m/mQ)| = 4.6 at m = 10" 



-18 



6.2. Comparison with existing studies 



We illustrate how our work fits into existing studies of fragmentation cascades with a parameter 
space plot in |U — /3 coordinates. Figure [5] shows the domains of applicability in the /i — /3 plane for 
the Dohnanyi and OBG solutions in relation to our analytic solutions for the monodisperse case. 
Each of our solutions is labeled by its corresponding formula number, and it is evident that our 
investigation covers the remainder of the jU — /3 plane. 




Fig. 5. — Pa r amet er space plot in the /U — /3 plane. The case considered bv iDohnanvil (119691): 
Tanaka et al.l (j 19961 ) is at the point /U = 1, /3 = 1, and the case considered bv lO'Brien k, Greenbere 
(|2003l ) lies on the line /i = 1 (solid line). Our solutions cover the rest of phase space and are 



labeled with references to corresponding equations in this work. Thus, solutions (IB25P and ()B23p 
(monodisperse, m_B(m) = Bm) lie on rays < 1, (3 = 1 and /i > l,/3 = 1, correspondingly. 
Solutions (IB26P and (IB27P (monodisperse, mB{m) = Bm^) are valid in the half planes /3 > 1 
(white) and /3 < 1 (gray) respective l y. Fo r these solutions ip =const along the line = 1 in 
agreement with lO'Brien &: Greenbergj (j2003l ). 



We next discuss why previous authors have not seen non-power law behavior. The reason is 
that they have all assumed mo,_B oc nit, and we have shown in §3.31 that non-power law behavior 
only results when mn,R is not prop orti onal to rrit. In some studies, the assumption mo,B oc rrit was 
explici t such as inlDohnanyil (119691) andlO'Brien &: Greeiiberd (120031 ) who both assumed mn = Cnit, 



and in IPetit fc Farinellal (1199 



O'Brien Sz Greenbers 



, and Ide Elfa &: Bruninil (|2007l ) who 
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all ass umecj^ mo^B = Trnt/2. In other cases, such as Tanaka et al. ( 1996 ) and Kobavashi &: Tanaka 
([20101) the scaling rriQ B on mt was implicit in assumptions about the form of fej ( H3.1\^3^ . 



We now return to our argument from ^2.31 that ruQ ^ should have the form (1371). T his con- 
clusion was based on an extension of the experimental results of iFuiiwara et alJ (119771^ beyond 
the strength-dominated regime. We mention t hat a number of authors (jPetit &: Farinellal . 119931 : 
O'Brien fc: Greenberei . l2005l : Ide Eha &: Bruninil . 120071 ) have considered a different extension of those 
results, and instead of our Eq. ([3T]) . these authors used 



mo{mt,mp) 



oc 



mt 



Ecou{mt,mp)/2 
Qs{'mt)mt 



(87) 



Here, Qs is the energ y per unit mass required to sh atter an object, but not necessarily to disperse its 
fragments to infinity (jO'Brien &: Greenbera . |2005| ) . If we assume for simplicity that oc Qs, then 
from Eq. (j87p and Eq. ()20p . we have tuq^b oc m^, for any functional fo rm of mR. The reason it is 
possible to derive two different forms for mo(m(, rrip) from the results of lFuiiwara et al.l (|l977l ) (Eq. 
()3ip and Eq. ()87p ) is because their experiments were performed over a small range of target masses 
using a constant projectile mass. The question of how to properly exte nd their results over a. larger 



200S 



mass range can best be settled by more experiments and simulations ([Stewart &: Leinhardt 
Benz &: Asphaugj . Il999l : iBenz et all . ll994l ). which can decisively answer how mo,s varies with mass. 
However, we point out that unless mo,B is exactly proportional to m^, non-power law behavior 
will result. As we show in the next section, these deviations from power law behavior can be 
observationally significant when extrapolating over many orders of magnitude in mass. 



6.3. Applications 

Our results clearly demonstrate that one should be somewhat cautious when adopting a pure 
power law approximation to describe the properties of fragmentation cascades. Even though the 
non-power law corrections computed in this work scale very weakly with object mass (as the square 
root of the logarithm of the mass (Eq. ([84p )). one has to keep in mind that in astrophysical systems 
collisional cascades span ~ 30 orders of magnitude in mass. Thus, even a weak deviation from a 
power law can become important, such as when inferring the disk mass ( dominated by the largest 



bodies) from its infrared luminosity (dominated by the smallest bodies) ([Wvattl . l2008l ) . 



Just for illustration, let us consider a population of Rinj = 10 km objects which get ground 
down to Rrm =1 Aiiii size particles by collisions. Infrared observations give us some idea of the 
mass in small particles, thus fixing the normalization of the mass spectrum at its low-mass end, 
and we want to infer from these data the total mass in large bodies feeding this collisional cascade. 
Connecting the mass contained at the low and high mass ends of the spectrum by a simple power 



^Note that mo,s — mt/2 does not follow from mrem{mB{'mt),mt) = mt/2, because the remnant does not belong 
to the distribution of ejecta O. 
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law leads to an error caused by the neglect of the non-power law effects. We can estimate this error 
6 by using Eq. (j84p and taking the ratio of ip at the high and low mass ends. Assuming some values 
of n and C in Eq. ([37|) that are "averaged" over the whole cascade (in practice these parameters 



will change several times between Rr 
objects), we have 



and Rinj because of variations in the internal properties of 



/ in(i;,^/i;*) 

^r^{Rinj/Ro) 



rrn 

\n{RQ/ Rinj) 



where Rq is the radius of the object with mass mg defined in §5.11 Assuming for illustration 
that on "average" /x ~ 1.1 and that at the high mass end the largest fragments produced in 
collisions have mass equal to 0.3 of the target mass (mo,B(?7imj) = 0.3mj„j) one finds ln{RQ/ Rinj) = 
ln(mo,_B(rraj„j)/mj„j)/3(;U — 1) f« 4 and 6^2 — 3. Thus, in this particular exercise the neglect of 
non-power law effects leads to an underestimate of the mass in large bodies by a factor of several. 
This also implies that the total disk mass, which is dominated by the mass in large bodies, is 
underestimated. 

An underestimate of the disk mass also leads to an underestimate of the disk lifetime, M^j^^ /M^j^^ . 
This occurs because if the disk is in steady-state, then F{m) is independent of m, which means it 
i s possible to infer M^isk from infrared observations alone, without extrapolation to large masses 
(jWvattl . I2OO8I ). Supposing that we had correctly inferred M^isk from observations, but had failed 
to apply the non-power law correction, and hence underestimated the disk mass, then we would 
also have underestimated the disk lifetime. 

Based on the above discussion, breaking the assumption tuq^b oc mt affects the calculation of 
disk properties from observations. Conversely, observations of disks can be used to constrain the 
model parameters (i.e. /i and C if mo,_B is given by Eq. (j37p ). if e.g. the inferred disk mass is found 
to be unreasonable for some parameter range. However, as mentioned in ^6.21 direct application of 
our theoretical results to the observed mass spectrum of objects is complicated by the multitude 
of additional factors playing an i mpor t ant role in real astrophysical systems. Nevertheless, modern 



calculations (de Eha &: Brunini 



2007 



O'Brien &: Greenbere 



2005 ) of the collisi o nal e volution in 



the asteroid belt and of debris disks (iThebault &: Augereaul . 120071 : iKrivov et al.l . l2008l ) aim for a 
precision of tens of percent or less over a broad range of masses. At this level of accuracy, the 
non-power law effects considered in this work would play a significant role and should be taken into 
account. 



7. Summary 

We have shown that unless raQ^B{fnt) on mt, where mo,B is the mass of the largest fragment 
produced in a collision with just enough energy to disperse half of the target's mass, mt, to infinity, 
a steady-state power law solution for the mass distribution, n{m), is not possible. The non-power 
law behavior is weak, however, and the solution for n[m) becomes the product of a power law and 
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a much more slowly varying function of the mass: n{m) = m "'(p{m), \d\inp / din m\ <^ 1. This 
slowly varyi ng function is equal t o a constant when mnR(rn t ) oc nit, and the fact that previous 



researchers (Kobavashi Sz Tanaka . 201C)I : Tanaka et al. . 199^: Dohnanvi . 1969 : Petit &: Farinella . 



I993I : lo'Brien k GreenbertJ . bood . boosl : Iwilliams fc Wetherilll . Il994l : Ide Elfa k Bruninil . boO?! ) have 



assumed just such a dependence of mo,_B(mt), explains why this kind of non-power law behavior 
was not observed earlier. 

When mQ^B{fnt) is not proportional to mt, n{m) deviates smoothly away from a pure power 
law, with the deviation only becoming significant when considering many orders of magnitude 
in mass. This is quite different from the wavy non-power law behavior resulting from either 
a low er cutoff to the mass distribution due to the eject i on of small particles by radiat ion pres- 



sure (IThebault &: Augereaul . 



2007 



Campo Bagatin et al 



1994 



200c 



Durda fc Dermottl. 119971). a tran- 
sition from a collisio nally-evolved to a primordial size distribution ( Fraseil . 



Pan &: Sari. 



2005- 



Kenyon fc Bromlevl . 12004). or a break in the power law index of the strength law (jo'Brien &: Greenber^ . 



2003 



2OO5I ). The non-power law behavior we describe in this work is significant when extrapolating 



over many orders of magnitude in mas s, such as wh en inferring the number of large bodies in a 



system based on infrared observations (jWvattl . 120081 ). For instance, assuming niQ^Bi^t 



oc m 



1.1 



deviation of only 10% in the power law index from the usual assumption of w.o,_B(mt) oc mt results 
in a factor of ~ 2 — 3 correction when inferring the number of 10 km bodies from observations of 
dust. 

We have quantified precisely the effect of the non-power law behavior on the mass distribution 
by obtaining analytical solutions for ip{m) in the case of a power law fragment mass distribution 
with mB{mt) = mt, and a monodisperse fragment mass distribution (all fragments the same size) 
with mB{mt) = Bm^ . We have also provided a general framework for solving the mB{mt) = Bm^ 
case with an arbitrary fragment mass distribution. In all cases considered, our analytical solutions 
were confirmed numerically, and we noticed that the simple expression (j84p captures the essence of 
the non-power law behavior for a wide range of parameters in our model. 

In the course of our investigation, we have also found an entirely different type of non-power 
law behavior. Namely, we have discovered that fragmentation cascades can support wavy, steady- 
state solutions, even when there is no upper or lower mass cutoff, and the strength law is given by 
a pure power law. Our results were derived for the monodisperse case, but such waves may also 
be able to persist for more realistic fragment mass distributions. In astrophysical systems, these 
kinds of waves could be triggered in stochastic collisions between large planetesimals that generate 
enough collisional debris to significantly alter n(m). 

We are grateful to Jeremy Goodman for useful discussions. The financial support for this work 
is provided by the Sloan Foundation and NASA via grant NNX08AH87G. 
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A. Convergence of the Mass Flux Integral 

We discuss here the convergence of the mass flux integrals in ^3.113.21 If the correct value of 
a is substituted into Eq. (j42p or Eq. ()43p . then the mass flux takes the form 

1 oo 

F{m) (X J dxx^^ J dyy^'^ f {x , y) . (Al) 



It is helpful to discuss under what circumstances the integrals in this expression converge. 

The integral over y converges at its upper limit if a > 1, because /(x, y) < 1. There is actually 
already a more stringent condition of a > 5/3, which comes from the requirement that the cross- 
section be dominated by the smallest particles ( §2.ip . so the requirement a > 1 is automatically 
fulfilled. 

We now consider convergence of the integral over y at its lower limit. Generally speaking, the 
quantity f{x,y)y~'^ drops off for y <C 1, the threshold for catastrophic breaking, which leads to 
convergence. In fact, in a model without erosion, f{x,y) = for y < 1, so there is a sharp cutoff 
at y = 1. In reality, erosion will smooth this cutoff, but as long as f{x,y) falls off faster than y"^^ 
for y <C 1, then the integral will converge at its lower limit. We point out that if the integral over y 
converges, then it follows from the above arguments that f{x,y)y~°^ is peaked at ypeak = k. In the 
absence of erosion, it i s clea r that A; ~ 1, but when erosion is considered, we typically have A; <C 1 



(|Kobavashi Tanakal . l20ld ') 



We next consider the convergence of the integral over x, and we find it helpful to define 

/>oo 

f{x)^ / dyy--/(x,y) (A2) 



As long as f{x) is bounded on x G [0,1], then the integral over x in Eq. ()Aip converges at the 
upper limit of integration. At the lower limit of integration, the function x~^ diverges, but only 
logarithmically, so if /(x) — > for x <C 1, then we would typically expect convergence at the lower 
limit. From Eq. ([7]) and Eq. (f22]l . we have /(x, y) — )■ for x <C 1, so we do indeed expect /(x) 
in the same limit. 



This preprint was prepared with the AAS lATpjX macros v5.2. 
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B. Details of the calculation of ip{m). 

We start by developing a general method for calculating 99(771) for different forms of mB(m) 
from the functional Eq. (j57p with known Q(m). 



B.l. Case mB{m) = Bm, /3 = 1 

Here, we will first assume following Dohnanyi (1969) and Tanaka et al. (1996) that msim) = 
Bm. Then, we need to solve the functional equation 

(p{m)(f{Bm) = Q{m). (Bl) 

Taking the logarithm of both sides of this equation we get 

lnip{m) + Imp (Bm) = In @{m), (B2) 

and upon introducing the new independent variable v = lnm, the constant b = InB, and the new 
function (pi{s) = lnip{e^) one gets 

ipi{v) + ipi{v + b) = Ri{v), Ri{v)=lne{e''). (B3) 

For some applications, it is more appropriate to study an equivalent equation 

ipi{u-b) + ipi{u) = Riiu-b). (B4) 



One can check by direct substitution that the formal solution of Eq. (jBSp up to an additive 
constant 921 (—00) is given by 

00 

Mv) + <^i(-oo) = Y,i-'^)''Mv + kb), (B5) 

fc=0 

where we have used the fact that 6 < 0, since B < 1. It follows then, that (p is given up to an 
overall normalization constant f{0) as 



(p{m) 

exp 



X](-l)'=lne(e'' 



fr 9(me^^'') 

ii e(me(2^'+i)^)- ^^^^ 



^In m+kb\ 
— ij mci^r 

.fc=0 

Analogously, the solution of Eq. ()B4p up to an additive constant and the corresponding solution 
for (/9 up to a normalization are given by 



^liv) + 991(00) = ^(-l)'^fli(^; -ik + 1)6), (B7) 



fc=0 



(^(00) e (me-(2fc+2)6) ■ ^ ' 

In TO. 3. II we provide an example of how to discriminate between using Eq. (|B6P versus Eq. (jBSP 
to calculate (p. 
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B.2. Case mBini ) = Bm^, /3 ^ 1 

Whenever /3 7^ 1 we need to solve the functional Eq. ip{m)(p (^Bm^^ = G(m), which is easily 
converted to 

lnip{m) +lnip (^Bm'^^ =lne{m). (B9) 

Invoking the mass scale m*^ = B^^^^~^^ at which breaking stops (i.e. mBim^) = m^), and 
defining the new independent variable u = In | ln(m/mg)|, — 00 < u < 00 and the new function 
(P2{s) = lnip (m^je'^'') one obtains 

ip2{u) + Mu + a) = R2{u), ii;2(n) =lnG(mj5e^") , a = ln/3. (BIO) 

This equation can also be converted to an equivalent form useful in some applications: 

(P2{u - a) + (p2{u) = R2{u - a). (Bll) 

Analogous to the previous case, we can write down formal solutions of Eq. (jBlOp and Eq. (|Blip 
up to additive constants, and the corresponding solutions for ip up to normalization constants. For 
/3 > 1, we have a > 0, which yields for Eq. (|B10p 

00 

Mu) + V'2(oo) = ^(-l)'=i22(n + ka), (B12) 

^ ^ 00 e^exp(e^^-ln^)) ^ ^^^^^ 
fi^) H © (m^j exp (^e('^k+i)a In^^^' 



and for Eq. (IBTTI) 



Mu) + M-00) = ^(-l)^i?2(n -{k + l)a), (B14) 



fc=0 



^(^) _ - G(m^^exp(e-(2'=+i)'^ln^ 



h ^ H\- (^1^) 



In a similar fashion, we can obtain the solutions when /3 < 1, in which case a < 0. We have 
for Eq. (|BT0]1 

CO 

ip2{u) + f2{-oo) = ^(-l)*'i?2(u + ka), (B16) 

fc=0 

^jm) - e(n.|,exp(e2^-ln^)) 
'^("^s) fcio e (m^ exp (^e{2fc+i)'^ In ^ ^ ^ ' 
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and for Eq. (|BTT]) 



oo 



V2{u) + 9?2(oo) = Y,i-^)^R2{u -{k + l)a), (B18) 

A:=0 



44 = n (bid) 



This completes the description of the general mathematical formalism needed for finding (^(m) 
given G(m) and given ms(m) in power law form. We now obtain explicit expressions for ip{m) for 
a monodisperse fragment mass distribution. 



B.3. Application to the monodisperse case 
B.3.1. Case mB{m) = Bm^ , p = I 
In the monodisperse case, we can take 

Q^"") = w / ^ ( = Tx \r TT^' (^20) 

ln(m/mo,B(mjj (1— ^jinm — InG 

where we will again assume ra^^^B^rn) = Cm^. Considering the case /3 = 1, the two relevant 
equations for obtaining (/?(m) are Eq. ()B6P and Eq. (jBSp . The 95(0) term in the denominator on 
the left hand side of Eq. (|B6P means that this equation is only applicable when > 1, since for 
^ < 1 there is a value rriQ below which mo^Bi^) > "i, and the solution is unphysical below this 
point. Similarly, the (p{oo) term in the denominator on the left hand side of Eq. (|B8p means that 
the solution only works for < 1, because for > 1, mo,B(m) > m above ttiq, and the solution is 
again unphysical. 

Treating first the // > 1 case, if we substitute Eq. ()B20p into Eq. ()B6p . we have 

00 In f + (2k + 1)(1 - n)b 

^ = n ^'"r ^ N —■ (B21) 

Unfortunately, this expression does not converge, implying that (/'(O) = 0- To understand this 
behavior, we can refer back to the analytic solution for ip in the monodisperse case with mBim) = m 
(Eq. ([El])). There, we had ip{m) = l/y^ln{m/mo^B{rn)). For > 1, this expression does indeed 
yield (/?(0) = 0. 

Since 99(771) is only defined up to a normalization, we can remedy the situation by working with 
the ratio of (p{m) at two points rather than with ip(m) itself, which gives a convergent expression. 
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If we define ai = ln(mi/mo,B(mi)), a2 = ln(m2/?TT-o, 5(771-2)), and c = (1 — /u)6, then it follows from 
expression ()B2ip that 

^{^1) ^ fx (ai + (2fc + l)6)(fl2 + 2A:fe) 

'/'("^s) l}^{ai + 2kb){a2 + {2k + l)by ^ ' 

which simplifies to 

virn,) ^ «2r(i + §^)r(i + f) 
v^M air(i + §|)r(i + f)- 

The correctness of this solution can be verified directly by substituting ip back into Eq. ([57|) . and 
using ([B20]) for e(m). 

Treating next the ji <1 case, if we substitute Eq. (jB20p into Eq. (jBSp . we have 



,^(^^ °o In — - (2fc + 2)(1 - p)h 

44 = IT y'""'"^"^^^ ^ —. (B24) 



Again, this equation does not converge, but the ratio of ip at two points does, and using the same 
definitions for oi, 02, and c as before, we have 

^ ^ r(i-f)rQ-^) 

r(i-|i)r(i-g)- ^ ^ 



5.5.5. Case mB(m) = Bml^ , (3^1 

We now treat the case when the mass of the smallest projectile that can catastrophically 
shatter a target is not proportional to the mass of the target itself, so that /3 7^ 1. 

We first consider the case /3 > 1, so we have the choice of using Eq. (jB13p or Eq. (jB15p . and we 
limit ourselves to the situation when mo,s(m.^) < m^. We expect such a situation to be realistic, 
since the maximum fragments created by the disruption of bodies with mass close to m^, should 
still be smaller than m^. Without this assumption, ip{m*^) would be unphysical, since we would 
have m'0,B(?n,^) > m^, which is impossible in a fragmentation cascade. Given this assumption, we 
can use Eq. (jBlSP to solve for 99 up to a normalization. 

Substituting Eq. (lB20]l into Eq. (IBT5]1 . we have 

,n(m\ lnf „ +ln£> 

K-u \^Dmo,fl(m) J ^ 

where we have defined the constant D = C^^B^^^^^^^^"^^ = (m^/mQ)^~^. In this case, the product 
does converge since ip{m*Q) is nonzero (it is also finite), and we do not have to take the ratio of two 
points as we did for the [3 = case in §B.3.1[ 
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Next, we consider the case /? < 1, and now we have the choice of using Eq. ()B17p or Eq. 
()B19p . We again make the assumption that mo^Bi'fn^) < w-^i which case we can use Eq. ()B17p 
to obtain 



ip[m) 



B) 



CO In 

n 

k=0 



-Dmo,s(m) 



+lnL» 



In 



Dmo,s(m) 



(B27) 



Again, there are no problems with convergence. Note, that if we set / x = 1 then m/mn,g(m) = cons t 
and Eq. ([B26]l and Eq. (|B27l) yield ip{m) =const in agreement with lO'Brien &: Greenbergi (j2003l ). 



C. Fragmentation Code 

We describe here the numerical algorithm we use to study fragmentation. We differ from the 
main text here in that our algorithm evolves the differential number density of particles n(r) per 
unit radius, rather than particle mass, but it is straightforward to convert between n(m) and n(r). 

The evolution equation for n(r) can be written as the sum of a source term (denoted by a plus 
sign) and a sink term (denoted by a minus sign): 

The sink term is simply given by the number of catastrophic collisions that particles with radius r 
are undergoing per unit time. Dropping the dependence on time for brevity (everything is evaluated 
at time t) and ignoring erosion, we can write 

— — (r) = —7rvn{r) / dr'n{r'){r + r') , (C2) 

dt Jrgir) 



where r^ir) is the minimum particle size that can fragment a particle of radius r, rBif) is defined 
analogously to fhB{m), and we have assumed that the particle cross-section is the geometric cross- 
section, and that the impact velocity is a constant §2.11 



For the source term, we will first state the equation and then analyze it piece by piece: 



9n+(r) 
dt 



fB(ri) 



TTV I drin{ri) / dr2n{r2)h{r,ri,r2){ri + ■ (C3) 

rB(ri) 



Here, we have used ri to denote a target and r2 a projectile with the possibility that r2 > ri. 
Now, TTvn{ri)n{r2){ri +r2)'^dr2dr\ gives the number of collisions (per unit volume per unit time) 
between targets in the range ri to ri -|- dri and projectiles in the range r2 to r2 + dr2- We define 
h{r,ri,r2)dr to be the number of fragments in the size range r to r + dr from destruction of the 
target particle only. Then h{r,ri,r2)n{ri)n{r2)TZ{ri,r2)dr2dridr gives the particle flux into the 
range r to r + dr from targets in the range ri to ri + dri that have been shattered by projectiles in 
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the range r2 to r2 + dr2. Then, we simply integrate this expression over all values of ri > r and all 
values of r2, which yield a catastrophic fragmentation event (i.e. an event in which both the target 
and the projectile are destroyed). 

A useful simplification of Eq. (IC3P can be obtained if we assume that the distribution of 
fragments is independent of the projectile size so that h{r,ri,r2) h{r,ri). In this case, which is 
adopted for the numerical calculations in this work (©, equation (jCSP becomes 



This is a useful form of the source equation when it comes to computations, because it reduces a 
double integral to a single integral, once the sink term has been calculated. 

Given a starting distribution for n(r) and a form for h{r, r') (e.g. monodisperse, power law, 
etc.), it is possible to evolve the distribution forward in time using discretized versions of Eq. 
()C2p and Eq. (jCSh in logr space. The integrals in these equations are performed using standard 
numerical integration techniques, such as the trapezoid or Simpson's rule. This results in an efficient 
order 0{N) method, where A'^ is the number of radius bins which are equally spaced in logr. The 
scaling of the method is important for us, since to resolve the non-power law behavior we use up 
to ~ 3000 bins per decade in r, yielding a total of ~ 10^ bins. 

In order to advance the distribution in time, a timestep must be specified. A good criterion 
is to set it to a fixed fraction of the shortest collision time in the simulation. This ensures that 
the particle size distribution can never become negative. Another problem is dealing with the 
boundaries of the simulation. The upper boundary is generally not problematic since the collision 
time there is long compared to the rest of the simulation, and the maximum particle size can be 
chosen to be as large as necessary for there to be negligible variation in the particle mass distribution 
at the high mass end. The lower boundary, however, can be a problem ii rsir) <C r there. In this 
case, extrapolation of the particle mass distribution to lower masses is necessary, and care must be 
taken to ensure that the simulation is stable. 




(C4) 



and comparing with Eq. IC2I we see that we can write 




(C5) 



